Inner edge of neutron-star crust with SLy 
effective nucleon-nucleon interactions 

F. Douchin^, P. Haensel^-^'^ 

^Centre de Recherche Astronomique de Lyon, 
Ecole Normale Superieure de Lyon, 46, allee d'ltalie, 
69364 Lyon, France 
"^N. Copernicus Astronomical Center, 
Polish Academy of Sciences, Bartycka 18, 
00-716 Warszawa, Poland 
^ Departement dAstrophysique Relativiste et de Cosmologie, 
UMR 8629 du CNRS, Ohservatoire de Pans, 
F-92195 Meudon Cedex, France 

February 1, 2008 



to he published in Physics Letters B 
Abstract 

The boundary between the sohd crust, assumed to be in complete thermodynamic 
equihbrium (cold catalyzed matter), and the liquid core of a neutron star is studied using 
Skyrme SLy effective N-N interactions. An approximate value of the density at the inner 
edge of the crust is obtained from the threshold for the instability of homogeneous npe 
matter with respect to small periodic density perturbations. Calculations of the bottom 
layer of the crust are performed within the Compressible Liquid Drop Model. Spherical 
nuclei are energetically preferred over exotic ones (cylinders, plates, tubes, bubbles), 
down to the inner edge of the crust, found at Pcdgc = 0.08 fm^'^. 
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The solid crust of a neutron star plays an important role in neutron star evolution 
and dynamics. It insulates thermally neutron star surface from its hot liquid interior, and 
therefore plays an essential role in neutron star cooling. It can build-up stresses, leading 
to instabilities responsible for the glitches in pulsar timing. Moreover, it can support non- 
axial deformations which, combined with rapid rotation, could make neutron star a source 
of continuous gravitational radiation. 

The mass of the neutron-star crust depends sensitively on the density of its inner edge 
and on the neutron-star matter equation of state. Recent calculations of this quantity 
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[17, 11, 18, 8] give Pedge — 0.08 — 0.10 fm~^, substantially lower than the estimate close to 
the normal nuclear density po = 0.16 fm~^, derived in the classical paper of Baym et al. [1]. 
Some theoretical calculations predict exotic distributions of nuclear matter (rods, plates, 
etc.) within the bottom layers of the crust [17, 11, 9]. However, the presence of exotic nuclei 
depends on the effective nucleon-nucleon (N-N) force used; for some forces, spherical nuclei 
are present down to the inner edge of the crust [17, 8]. Presence of exotic nuclei would lead 
to different transport and elastic properties of the crust, compared to the standard case of 
a Coulomb crystal formed by spherical nuclei [17, 13]. This, in turn, would have important 
consequences for dynamics and evolution of neutron star. 

In order to make calculations of the structure of the bottom layers of the neutron- 
star crust feasible, one has to use an effective N-N interaction. Frequently used in as- 
trophysical applications SkM* force was constructed to provide a consistent description of 
isovector properties of nuclei (giant dipole resonances) [5], and seems thus to be partic- 
ularly suitable for the description of neutron rich nuclei. However, this as well as many 
other existing effective N-N interactions were fitted to the properties of laboratory atomic 
nuclei, with (A'^ — Z)/A < 0.3, while in the bottom layers of neutron-star crust one ex- 
pects {pn — Pp) /p > 0.8. In view of this, application of these effective N-N interactions to 
the bottom layers of a neutron-star crust involves a rather risky extrapolation to strongly 
asymmetric nucleon matter. In order to remove a part of this uncertainty, modifications 
of effective N-N forces, to make them consistent with microscopic calculations of neutron 
matter, have been applied. Such a procedure was used in seventies to obtain the Ski' force 
[4], via a rather ad hoc modification of the Ski force constructed originally by Vautherin 
and Brink [3] to describe terrestrial nuclei. In this way, it became consistent with energy 
per nucleon of neutron matter calculated by Siemens and Pandharipande [15]. Later, gen- 
eralized types of the Skyrme interaction, FPS [12] and FPS21 [18], with larger number of 
fitted parameters and more general density dependence, were derived by fitting the temper- 
ature and density dependent energies per baryon of nuclear and neutron matter obtained 
in microscopic calculations of Friedman and Pandharipande [10]. 

Recently, a new set of the Skyrme-typc effective N-N interaction has been derived, 
based on an approach which may be more appropriate, as far as the applications to a very 
neutron rich matter are concerned [6, 7]. Relevant additional experimental items concerning 
neutron rich nuclei (isovector effective masses) , constraints of spin stability, and requirement 
of consistency with the UV14H-VIII equation of state (EOS) of dense neutron matter of 
Wiringa et al. [20] foi po < p < 1.5 fm~^, were combined with general procedure of fitting 
the properties of doubly magic nuclei. This procedure led to a set of the SLy (Skyrme Lyon) 
models, which - due to the emphasis put on their neutron-excess dependence - seem to be 
particularly suitable for the calculations of the properties of neutron-star crust. 

In the present Letter, we calculate the structure of the bottom layer of the neutron- 
star crust, and the density at its inner edge, Pedge; using the SLy models of effective N-N 
interaction, and compare our results with those obtained using older Skyrme-type forces, 
SkM* and Ski'. The parameters of the SLy forces used in the present calculations, together 
with those of the SkM* and Ski' models, are given in Table 1. The SLy4 interaction is a 
basic SLy force. The SLy7 model has been obtained following the most ambitious fitting 
procedure, in which spin-gradient terms and center of mass correction term were both 
included in the Skyrme energy functional [7] . 
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Recent calculations show that neutron-star matter (both liquid and solid) close to the 
crust-liquid core interface contains only a few percent of protons. Therefore, an effective 
N-N interaction used for the calculation of the crust-liquid transition should, in the first 
place, yield a realistic description of pure neutron matter at subnuclcar density. As wc do 
not have a direct access to the EOS of pure neutron matter at subnuclear densities, results 
of precise numerical calculations of the ground state at p < po, carried out using modern 
many-body theories and the best bare N-N hamiltonians, are to be used as an ersatz of such 
experimental data [18]. 

The SLy forces have been constructed as to be consistent with the UV14-I-UVII model 
of Wiringa et al. [20] of neutron matter above po [6, 7]. It is therefore of interest to see 
how well these effective N-N forces reproduce the UV14-I-UVII model of neutron matter at 
subnuclear densities. 

In Fig. 1 we have plotted the energy per neutron, e, versus neutron density, p < po, 
for the SLy4 and SLy7 effective N-N interactions. The filled squares correspond to realistic 
UV14-I-UVII model of neutron matter, and were taken from Table III of Wiringa et al. [20]. 

The agreement of the SLy4, SLy7 curves with the neutron matter EOS of the UV14+VII 
model of Wiringa et al. [20] is very good. This is not the case for older Skyrme forces. Ski' 
and SkM*. 

In what follows, a mixture of neutrons, protons, and electrons will be referred to as npe 
matter. An electrically neutral npe matter in beta equilibrium corresponds to neutron-star 
matter at p ^ pg. At a given p, the ground state of a homogeneous npe matter corresponds 
to the minimum of the energy density E{pn, Pp, pe) = Eq, under the constraints of fixed 
baryon density and electric charge neutrality, pp + Pn = P and pe = Pp, respectively. This 
implies beta equilibrium between the matter constituents and ensures vanishing of the first 
variation of E due to small perturbations 5pj{r) (where j = n,p,e) of the equilibrium 
solution (under the constraints of constant total nucleon number, A = pV, and global 
charge neutrality within the volume V of the system). However, this does not guarantee 
the stability of the spatially homogeneous state of the npe matter, which requires that the 
second variation of E (quadratic in 6pj) be positive. 

Our expression for the energy functional of slightly inhomogeneous neutron-star matter 
has been calculated using the semi-classical Extended Thomas-Fermi (ETF) treatment of 
the kinetic and the spin-gradient terms in nucleon contribution to E [2]. Assuming that 
the spatial gradients are small, we keep only the quadratic gradient terms in the ETF 
expressions. Our approximation is justified by the fact that characteristic wavelengths of 
periodic perturbations turn out to be much larger than the internucleon distance. With 
these approximations, the change of the energy (per unit volume) implied by the density 
perturbations can be expressed, keeping only second order terms [1, 18], 

^'^'^Uwf^ F,,(q),5ft(q),5p,(q)* , (1) 

where we used the Fourier representation 

Spj{r) = J J^jsSpMW''^ ■ (2) 

The hermitian Fij (q) matrix determines the stability of the uniform state of equilibrium of 
the npe matter with respect to the spatially periodic perturbations of wavevector q. Due 
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to the isotropy of the homogeneous equihbrium state of the npe matter, F^j depends only 
on |q| = q. In the case of Skyrme-type forces, the matrix elements Fij can be calculated 
analytically, as explicit functions of the equilibrium densities and q, and are composed of 
compression (local), curvature (gradient), and Coulomb components [1, 18]. 

The condition for the F^j matrix to be positive definite is equivalent to the requirement 
that the determinant of the Fij matrix be positive [18]. At each density p, one has thus to 
check the condition det[Fij{q)] > 0. Let us start with some p, at which det[Fij{q)] > for 
any q. By decreasing p, we find eventually a wavenumber Q at which stability condition is 
violated for the first time; this happens at some density pQ. For p < pg the homogeneous 
state is no longer the ground state of the npe matter since it is unstable with respect to 
small periodic density modulations. 

We get pq = 0.079 fm~^ for both SLy4 and SLy7 forces, to be compared with 0.075 fm~^ 
and 0.101 fm~^ for the SkM* and Ski' forces, respectively. The instability at pQ signals 
a phase transition with a loss of translational symmetry of the npe matter, and appear- 
ance of nuclear structures. Calculations performed previously with Skyrme forces indicate 
that Pq ~ Pcdgc [18]. The differences between the values of pQ between four forces are 
correlated with differences in the density behavior of nucleon chemical potentials, whose 
density dervatives enter the stability matrix F^j [18]. In Fig.2 we plot the values of pn and 
Pp for all four Skyrme forces considered, versus nucleon density p, for the npe matter in 
beta equilibrium. In the relevant density region, p ~ ^po, one notices very good agreement 
between the values of pp. The relative differences in the values of and in particular in 
the slopes of the Pn{p) curves are much larger. The slope of Pn{p) for the SkM* force is the 
highest, and that for the Ski' force the lowest, at p ~ ^pQ. Notice that the slope of pn for 
the Ski' is quite similar to those obtained in [18] for the EPS and FPS21 forces (see Fig. 5 
of [18]). Qualitatively, the higher the slope of the lower the value of pq (c.f., [18]). 

Nuclear structures in bottom layers of neutron-star crust were described using the Com- 
pressible Liquid Drop Model (CLDM) of nuclei [1, 21, 22, 23]. Within the CLDM, one is 
able to separate bulk, surface, and Coulomb contributions to total energy density, E. Elec- 
trons are assumed to form an uniform Fermi gas, and yield the rest and kinetic energy 
contribution, denoted by E^. Within CLDM, total energy density of the inner neutron-star 
crust is given by 

E = £^N,bulk + -E^N,surf + -Ecoul + Eg . (3) 

Here, -En, bulk is the bulk contribution of nucleons, which does not depend on the shape of 
nuclear structures. However, both i?N,surf and E'coul do depend on the shape of nuclear 
structures, formed by denser nuclear matter and the less dense neutron gas (detailed de- 
scription of the calculation of £'N,surf and -E/coui for the SLy forces was presented in [23] ; our 
results for the nuclear surface and curvature properties for SLy forces will be described in de- 
tail elsewhere [24]). We restricted ourselves to three shapes of the nuclear matter - neutron 
gas interface: spherical, cylindrical, and plane. Consequently, we considered five types of 
nuclear structures: spheres of nuclear matter in neutron gas, cylinders of nuclear matter in 
neutron gas (rods), plane slabs of nuclear matter in neutron gas, cylindrical holes in nuclear 
matter filled by neutron gas (tubes), and spherical holes in nuclear matter filled by neutron 
gas (bubbles). In view of a significant neutron excess, the interface includes neutron skin, 
formed by neutrons adsorbed onto the nuclear matter surface. In view of a finite thickness 
of nuclear surface, the definition of its spatial location is a matter of convention. Here, we 



4 



defined it by the radius of the equivalent constant density proton distribution, Rp, which 
determines thus the radius of spheres, bubbles, cylinders and tubes, and the half-thickness 
of plane nuclear matter slabs. The neutron radius, was defined by the condition that 
it yields a squarcd-off neutron density ditribution with constant neutron densities, which 
are equal to the real ones far from the nuclear matter - neutron gas interface, reproduces 
actual numbers of neutrons. The thickness of the neutron skin is then defined as Rn — Rp, 
Rp being defined in a similar manner as Rn . The nuclear surface energy term, -ENjSurf, 
gives the contribution of the interface between neutron gas and nuclear matter; it includes 
contribution of neutron skin [21, 16, 22]. In the case of spherical and cylindrical interface, 
£^N,surf includes curvature correction; the curvature correction vanishes for slabs. 

In order to calculate £^Coui5 we used the Wigner-Seitz approximation. In the case of 
spheres, bubbles, rods, and tubes, Wigner-Seitz cells were approximated by spheres, and 
cylinders, of radius -Rceii- In the case of slabs, Wigner-Seitz cells were bounded by planes, 
with Rceii being defined as the half-distance between plane boundaries of the cell. At given 
average nucleon density, p, and for an assumed shape of nuclear structures, the energy 
density was minimized with respect to thermodynamic variables, under the condition of an 
average charge neutrality [21, 22, 23]. Our results for E{p) for the SLy4 force, together 
with the energy density of the uniform npe matter, are displayed in Fig. 3. For the sake 
of convenience, energy densities has been shown after subtracting the value of the energy 
density of the bulk nuclear matter - neutron gas system (charge neutral and in beta equilib- 
rium), AE(p) = -E'(p) — -Ebuik(p)- Spherical nuclei turn out to be energetically preferred over 
other nuclear shapes, and also over homogeneous npe matter, down to pcrit = 0.077 fm~^. 
Therefore, within our set of possible nuclear shapes, the ground state of neutron-star crust 
contains spherical nuclei, down to its bottom edge. 

As was already noticed in [28, 29, 17], exotic nuclei will be present only if the filling 
fraction, u, is large enough. For a given neutron excess, introducing finite size (Coulomb 
and surface) effects increases the value of u. This increase of u is larger when the surface 
tension is larger, the finite-size effect being roughly propotional to it. Therefore, the larger 
the surface tension, the more likely the exotic nuclei are to appear. For all Skyrme forces 
we have studied so far, our results are in agreement with this qualitative rule. In particular, 
at the proton fraction ~ 5 — 8% the surface tension for the SkM* and SLy forces, for which 
exotic nuclei do not appear, is much smaller than for the Ski' model, for which exotic nuclei 
do appear. 

Actually, under conditions of thermodynamic equilibrium, transition from the crust to 
the uniform liquid takes place at a constant pressure, and is accompanied by a density 
jump (first order phase transition). Using Maxwell construction, we find that the edge of 
the crust has density Pcdgc = 0.076 fm^^, and coexists there with uniform npe matter of the 
density higher by Ap/pedge = 1-4%. Therefore, crust-liquid core transition is a very weak 
first-order phase transition; it takes place at Pedge = 0.34 MeV/fm^. 

Geometrical parameters characterizing the bottom layers of the crust are shown in Fig. 
4. Here, Rp is the proton radiiis corresponding to the equivalent squared proton distribution 
within nuclei, i^ceii is the radius of the Wigner-Seitz cell, and u is the fraction of volume 
occupied by nuclear matter (equal to that occupied by protons). The crust-liquid core 
transition takes place not because nuclei grow in size, but rather because they become 
closer and closer. The filling fraction u grows rapidly for p approaching Pedge; from only 
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5 % at p = 0.04 fm~^ to nearly 30 % at Pedge- Still, less than 30% of the volume is filled 
by nuclear matter at the crust-liquid core transition point. Our maximum filling factor 
is significantly larger than the limiting value of 0.2 for spherical nuclei phase, derived by 
Oyamatsu et al. [27] using general considerations involving surface and Coulomb energy. 
This is due to the curvature corrections, included in our formula for -EN sm-f and absent 
in the model of Oyamatsu et al. [27]. Curvature tension is particularly important near 
the crust-liquid transition point, where the surface tension is very small. If we neglect the 
curvature term in surface energy, we find maximum u for nuclei of about 0.2, consistent 
with Oyamatsu et al. [27]. Finally, let us mention that no proton drip occurs in the ground 
state of the crust. 

More detailed information on neutron-rich nuclei, present in the ground state of the 
bottom layers of neutron-star crust with p > ^po, can be found in Fig. 5. Number of 
nucleons in a nucleus. A, grows monotonically, and reaches about 600 at the edge of the 
crust. However, the number of protons changes rather weakly, from Z ~ 40 near neutron 
drip, to Z c^50 near the edge of the crust. Our results for Z of spherical nuclei are similar to 
those obtained in [11, 14], but are somewhat higher than those obtained using a relativistic 
mean-field model in [9]. An interesting quantity is the number of neutrons forming neutron 
skin, iVn^surf- For p > ^po, we find that A^n,surf decreases with increasing p, despite an 
increase of nuclear radius. This behavior is due to the fact, that neutron densities outside 
and inside nuclei become more an more alike, while thickness of neutron skin decreases, 
with increasing p, for p > |po- 

Near the bottom of the crust, spherical nuclei are very heavy, A ^ 300 — 500, and doubts 
concerning their stability with respect to deformation and fission arise. An approximate 
condition for fission in dense neutron-star matter is Rp/Rcdi > 1/2 [16]; this condition 
is fulfilled near the edge of the crust. Fig. 4. However, curvature terms in -EN.surf and 
-E'Coub which were not included in generalized Bohr- Wheeler condition of [16], as well as 
shell correction to E [26] , stabilize large spherical nuclei. Clearly, the problem of stability 
of spherical nuclei near the bottom of neutron-star crust deserves further investigation. In 
any case, however, forming rod-like nuclei in the ground state of neutron-star crust seems 
to be excluded by the energetical arguments. Fig. 3. 

In the CLDM approximation with SLy4 force, uniform npe matter is present in the 
ground state of the neutron-star interior at the density exceeding pnq = 0.077 fm~^. Com- 
paring the value of pnq with pg, we have an independent test of precision of our CLDM 
calculations. The relative difference between these two quantities is less than 3%, which is 
of the order-of-magnitude of the relative density jump in the crust - liquid core transition, 
obtained within the CLDM. For the SLy4 force, the crust - liquid core transition takes thus 
place at half of the normal nuclear density, and very similar value is obtained for the SLy7 
force. 

The structure of the ground state of neutron-star interior at subnuclear densities, deter- 
mined in this Letter for the SLy4 force, implies that the bottom layer of the neutron star 
crust is a three-dimensional Coulomb crystal, with spherical nuclei in the lattice sites, the 
edge of the crystal being located at Pedge — 0.08 fm^'^. 

We combined the SLy4 EOS of the inner crust with that of the liquid core, getting 
in this way a physically consistent, unified description of the neutron star interior beyond 
neutron drip. As for the outer crust, we have chosen the EOS of Haensel and Pichon [25], in 
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which maximum use of recent experimental information on neutron rich nuclei was made. 
For neutron star with canonical mass of 1.4 Mq, the crust contains AMcrust/-^ = 1-2% 
of the total stellar mass, and its relative contribution to the total moment of inertia is 
A/crust/-^ = 2.6%. The corresponding values for the FPS neutron star model arc some 
30% lower, despite larger value of Pedge(FPS) = 0.10 fm^-^. It turns out, however, that 
Pedge(FPS) ~ Pedge(SLy4), bccausc the FPS EOS near Pedge is significantly softer than the 
SLy4 one. The difference in the AMcmst/M for both models can thus be explained by 
the difference in the - 2GM/Rc^) factor appearing in an approximate expression for 
AMcrust [17]; the FPS neutron star of 1.4 Mq is more compact than the SLy4 one. Similar 
arguments can be used to explain the difference between the A/crust /-^ values for the SLy4 
and FPS forces (see [17] for the approximate formula to be used). However, the most 
important difference between the SLy4 and the FPS crusts is that in our case the whole 
crust, down to its bottom edge, is a Coulomb crystal, while in the FPS case the bottom 
layer, constituting about half of the crust mass, contains nonstandard nuclear shapes (rods, 
plates, tubes, bubbles), and therefore has elastic and transport properties very different 
from those of a three-dimensional Coulomb crystal [13]. Finally, presence of a Coulomb 
crystal is also relevant for the crust - liquid core boundary condition in rotating neutron 
stars. 
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Table 1 

Parameter values of the Skyrme forces 



force 


SLy4 


SLy7 


SkM* 


Ski' 


to (MeV im'^) 


-2488.91 


-2482.41 


-2645.0 


-1057.3 


ti (MeV fm^) 


486.82 


457.97 


410.0 


235.9 


t2 (MeV fm^) 


-546.39 


-419.85 


-135.0 


-100.00 


ta (MeV fm^+3'") 


13777.0 


13677.0 


15595.0 


14463.5 


(T 


1 

6 


1 

6 


1 

6 


1 


Xo 


0.834 


0.846 


0.09 


0.2885 


Xi 


-0.344 


-0.511 








X2 


-1.000 


-1.000 








X3 


1.354 


1.391 





0.2257 


Wo (MeV fm^) 


123.0 


126.0 


130.0 


120.0 
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Figure 1: 

Energy per neutron versus neutron density for neutron matter. Filled squares: results of 
Wiringa et al. [20] for the UV14+UVIII model. Solid line: the SLy4 model. Short dashes: 
the SLy7 model. For the sake of comparison, we show also results obtained for the SkM* 
(long dashes) and the Ski' (long dash-dot line) models. 
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Figure 2: 

Neutron and proton chemical potentials (with rest energy subtracted), Hn, iJ'p, versus total 
nuclcon density, p, for the SLy and SkM*, Ski' forces, respectively. The /Xp curve for the 
SLy7 force cannot be graphically distinguished from the SLy4 one. 
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Figure 3: 

Energy density of a given phase of inner-crust matter minus that of the bulk nuclear matter 
- neutron gas system, as a function of the average nucleon density, p. Thick solid line 
corresponds to the homogeneous npe matter. Calculations performed for the SLy4 force. 
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Figure 4: 

Radius of spherical Wigner-Seitz cell, Rcdi, the proton radius of spherical nuclei, Rp, and 
fraction of volume filled by nuclear matter, u (in percent), versus average nucleon density, 
p. Shaded area corresponds to the homogeneous npe matter. Calculations performed for 
the SLy4 force. 
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Figure 5: 

Mass number of spherical nuclei, A, and their proton number, Z, versus average nucleon 
density, p. Dotted line corresponds to number of nucleons after deducing neutrons belonging 
to neutron skin. Shaded area corresponds to the homogeneous npe matter. Calculations 
performed for the SLy4 force. 
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